library(tidyverse)
library(colorDF)
library(ggplot2)
library(ggbeeswarm)
library(tmod)
library(colorDF)
library(cowplot)
theme_set(theme_minimal())

A simple data set

Basic analysis

This is a data set that I used in my poster. It is available from GEO and contains RNA-Seq counts for healthy individuals and COVID-19 patients, described by Mick et al. (2020).

The data set has been already processed and downloaded; however, below you will find the code to replicate this procedure (if you want).

## this sometimes fails. The RDS of the object is provided
Sys.setenv(VROOM_CONNECTION_SIZE=8*131072)

## with getGEO, we can only get the phenoData
geo <- getGEO("GSE156063")[[1]]
saveRDS(geo, file="GSE156063.geo")
geo <- readRDS("GSE156063.geo")
covar <- as(phenoData(geo), "data.frame")

## remove columns with only one value
boring <- map_lgl(covar, ~ length(unique(.x)) == 1)
covar <- covar[ , !boring ] 

## clean up
covar <- covar %>% 
  dplyr::rename(gender = "gender:ch1") %>%
  mutate(disease = gsub(" .*", "", .data[["disease state:ch1"]])) %>%
  mutate(label = description) %>%
  mutate(group = disease) %>%
  arrange(description) %>%
  dplyr::select(all_of(c("title", "label", "gender", "disease", "group")))
## the counts must be downloaded from GEO separately.
if(!file.exists("GSE156063_swab_gene_counts.csv.gz")) {
  download.file("https://ftp.ncbi.nlm.nih.gov/geo/series/GSE156nnn/GSE156063/suppl/GSE156063_swab_gene_counts.csv.gz",
                "GSE156063_swab_gene_counts.csv.gz")
}

counts <- read_csv("GSE156063_swab_gene_counts.csv.gz")
.tmp <- counts[[1]]
counts <- as.matrix(counts[,-1])
rownames(counts) <- .tmp
#counts <- counts[ , covar$description ]
counts <- counts[ , covar$label ]
lcpm   <- edgeR::cpm(counts, log=TRUE)
#stopifnot(all(colnames(counts) == covar$description))
stopifnot(all(colnames(counts) == covar$label))

annot <- data.frame(ENSEMBL = rownames(counts))

.tmp <- mapIds(org.Hs.eg.db, annot$ENSEMBL, column=c("ENTREZID"), keytype="ENSEMBL")
annot$ENTREZID <- .tmp[ match(annot$ENSEMBL, names(.tmp)) ]

.tmp <- mapIds(org.Hs.eg.db, annot$ENSEMBL, column=c("SYMBOL"), keytype="ENSEMBL")
annot$SYMBOL <- .tmp[ match(annot$ENSEMBL, names(.tmp)) ]

.tmp <- mapIds(org.Hs.eg.db, annot$ENSEMBL, column=c("GENENAME"), keytype="ENSEMBL")
annot$GENENAME <- .tmp[ match(annot$ENSEMBL, names(.tmp)) ]
sel <- covar$group %in% c("no", "SC2")
counts <- counts[ , sel ]
covar  <- covar[ sel, ]
covar$group <- as.character(covar$group)

covar$group.disease <- paste0(covar$group, '_', covar$disease)
rownames(covar) <- covar$label
saveRDS(covar, file="covar.rds")
saveRDS(counts, file="counts.rds")
## DESeq2 calculations
## manual cache, since the operation takes a long time
sel <- !is.na(covar$group)

ds2 <- DESeqDataSetFromMatrix(counts[,sel], colData=covar[sel, ], design=~ disease)

ds2_file <- "ds2_cached.rds"
if(!file.exists(ds2_file)) {
  message("Running DESeq2")
  ds2 <- DESeq(ds2)
  saveRDS(ds2, file=ds2_file)
} else {
  message("Reading ds2 from manual cache")
  ds2 <- readRDS(ds2_file)
}

First, take a look at the structure of the data:

covar  <- readRDS("data/covar.rds") # covariates
summary_colorDF(covar)

The results table is stored under “data/disease_SC2_vs_no.tsv”. In data/annotation.tsv, there is the annotation of the genes in this data set.

dsres <- read_table("data/disease_SC2_vs_no.tsv")
annot <- read_table("data/annotation.tsv")
dsres <- merge(annot, dsres, by="ENSEMBL")

Basic CERNO analysis:

gl <- dsres$SYMBOL[ order(dsres$pvalue) ]
tres <- tmodCERNOtest(gl)
print(tres)
## # tmod report (class tmodReport) 8 x 107:
## # (Showing rows 1 - 5 out of 107)
##         │ID      │Title                  │cerno│N1   │AUC  │cES  │P.Value
##  DC.M3.4│ DC.M3.4│Interferon             │  425│   48│ 0.89│  4.4│<2e-16 
##  DC.M1.2│ DC.M1.2│Interferon             │  296│   24│ 0.99│  6.2│<2e-16 
## DC.M5.12│DC.M5.12│Interferon             │  356│   56│ 0.82│  3.2│<2e-16 
##  DC.M3.2│ DC.M3.2│Inflammation           │  528│  112│ 0.82│  2.4│<2e-16 
##   LI.M75│  LI.M75│antiviral IFN signature│  198│   21│ 0.91│  4.7│<2e-16 
##         │adj.P.Val
##  DC.M3.4│  3.9e-40
##  DC.M1.2│  5.6e-35
## DC.M5.12│  7.7e-25
##  DC.M3.2│  1.8e-24
##   LI.M75│  5.6e-20

Which we can quickly visualize:

ggPanelplot(tres)

We can also add information about whether the genes go up or down. Note that we need to put tres inside a list, so ggPanelplot can match tmod results with the sgenes object.

sgenes <- tmodDecideTests(dsres$SYMBOL, lfc = dsres$log2FoldChange,
                          pval=dsres$pvalue)
names(sgenes) <- "SC2_vs_no"
tres <- list("SC2_vs_no"=tres)
ggPanelplot(tres, sgenes=sgenes)

Using another gene sets

Another built-in gene set are the cell surface markers.

data(cell_signatures)
tres2 <- tmodCERNOtest(gl, mset=cell_signatures)

tres2 <- list("SC2_vs_no"=tres2)
sgenes2 <- tmodDecideTests(dsres$SYMBOL, lfc = dsres$log2FoldChange,
                          pval=dsres$pvalue, mset=cell_signatures)
names(sgenes2) <- "SC2_vs_no"
ggPanelplot(tres2, sgenes=sgenes2, mset=cell_signatures)

The mset parameter is used in most of the functions to use another gene sets.

Inspecting individual enrichments

Let us consider the gene set “LI.M68” – RIG-1 like receptor signalling. First, show the genes in this gene set:

data(tmod) # make the built-in gene sets visible
getModuleMembers("LI.M68", mset=tmod)
## $LI.M68
##  [1] "DDX58"  "IRF7"   "DHX58"  "TRIM25" "NFKBIA" "IFIH1"  "TNF"    "ISG15" 
##  [9] "CXCL10" "IL8"

Evidence plot – a ROC curve:

par(mfrow=c(2,2))
evidencePlot(gl, m = "DC.M1.2", gene.labels=TRUE)
evidencePlot(gl, m = "LI.M68", gene.labels=TRUE)
evidencePlot(gl, m = "LI.M165")
evidencePlot(gl, m = "LI.M4.0")

Eigengenes

The basic idea is as follos: run PCA only on a selected group of genes, use PC1 instead of average gene expression of that gene set (mind the signs!). Below, counts.rds contains the raw RNA-seq counts, which I first convert to log2 counts per million.

counts <- readRDS("data/counts.rds")

lcpm <- edgeR::cpm(counts, log=TRUE)
mm <- getModuleMembers("LI.M68", mset=tmod)[[1]]
sel <- annot$SYMBOL %in% mm
covar$eig <- prcomp(t(lcpm[ sel, ]), scale.=TRUE)$x[,1]

We can use eigengenes to nicely plot the results:

ggplot(covar, aes(x=disease, y=eig)) + geom_boxplot(outlier.shape = NA) +
  geom_beeswarm(cex=3)

(However, mind the sign!)

In tmod, eigengenes can be computed directly. Tmod makes sure that the eigengene correlates positively with the majority of the genes in the gene set and reverses the sign when necessary:

selm <- tres[[1]]$ID[1:4]
eig <- eigengene(lcpm, annot$SYMBOL)
eig <- as.data.frame(t(eig)) %>%
  select(all_of(selm)) %>%
  rownames_to_column("label") %>% 
  inner_join(covar, by="label")
plots <- map(selm, ~ ggplot(eig, aes_string(x="disease", y=.x)) +
             geom_boxplot(outlier.shape=NA) + 
             geom_beeswarm(cex=3) +
             ggtitle(tres[[1]]$Title[ match(.x, tres[[1]]$ID) ]))
plot_grid(plotlist=plots)

Combining gene set enrichments with PCA

Using a method like tmodCERNOtest, tmodUtest or fsgsea makes it possible to apply gene set enrichment to any single ranked list of genes.

(I have selected the components based on correlation with covar$disease)

pca <- prcomp(t(lcpm), scale.=TRUE)
pca_df <- pca$x %>% as.data.frame %>% 
  rownames_to_column("label") %>% merge(covar)
ggplot(pca_df, aes(x=PC5, y=PC6, color=disease)) + geom_point()

We can sort all the genes by their decreasing absolute loadings in the pca object.

pcares <- map(1:10, ~ {
                gl <- annot$SYMBOL[ order(-abs(pca$rotation[,.x])) ]
                tmodCERNOtest(gl)
             })
names(pcares) <- colnames(pca$rotation)[1:10]
sgenes <- tmodDecideTests(annot$SYMBOL, 
                          lfc=pca$rotation[, 1:10],
                          pval=NULL,
                          lfc.thr = 0)
ggPanelplot(pcares, sgenes=sgenes, filter_row_auc = .85)

ggplot(pca_df, aes(x=PC2, y=PC3, color=disease)) + geom_point()

Check the function tmod_pca for more details.

Using gene sets from msigdbr

Converting the data frame from the msigdbr package for use with tmod is straightforward:

library(msigdbr)
msig_df <- msigdbr("Homo sapiens")

goset <- msig_df %>%
  filter(gs_subcat == "GO:BP") %>%
  makeTmodFromDataFrame(module_col="gs_name",
                        title_col="gs_name",
                        feature_col="gene_symbol")
library(DT)
tmodCERNOtest(gl, mset=goset, qval = 1e-3) %>%
  filter(N1 > 5) %>%
  datatable() %>%
  formatSignif(c("adj.P.Val", "P.Value")) %>%
  formatRound(c("cES", "AUC"))

Upset plots

tres[[1]] %>% filter(N1 < 20) %>% pull(ID) %>% upset()

Overlaps between gene sets

foo <- tres[[1]] %>% 
  slice(1:20) %>%
  pull(ID) %>% 
  modCorPlot(stat = "overlap")
foo

“Leading edge” analysis

Which genes are driving the enrichment?

evidencePlot(gl, m = "LI.M165", mset=tmod)

evidencePlot(gl, m = "LI.M165", mset=tmod, style="gsea")

Get the leading edge:

lea <- tmodLEA(gl, "LI.M165")
lea
## $LI.M165
##  [1] "IFIT1"    "IFI27"    "CXCL11"   "RSAD2"    "SIGLEC1"  "PARP9"   
##  [7] "LAMP3"    "IFIT2"    "IFIT3"    "IFIH1"    "HERC5"    "CCL8"    
## [13] "SERPING1" "HESX1"    "CCL20"    "DSE"      "IFNGR2"   "RBM47"   
## [19] "MX2"      "BCL2A1"   "TNFAIP6"  "FAM20C"   "CSF2RB"   "DOCK4"   
## attr(,"N")
## [1] 35
## attr(,"LEA")
## [1] 2831
## attr(,"LEA.frac")
## [1] 0.6857143
## 
## attr(,"class")
## [1] "tmodLEA" "list"
clr <- ifelse(gl %in% lea[[1]], "red", "grey")
names(clr) <- gl
evidencePlot(gl, m = "LI.M165", mset=tmod,
               gene.labels=TRUE,
               gene.colors = clr)

And GSEA style:

clr <- ifelse(gl %in% lea[[1]], "red", "grey")
names(clr) <- gl
evidencePlot(gl, m = "LI.M165", mset=tmod,
               gene.labels=TRUE,
               gene.colors = clr,
               style="gsea")
abline(v=attr(lea[[1]], "LEA"), col="red", lwd=2)

The vaccination example

This example data set has been published by Weiner et al. (2019). RNA has been collected daily from individuals vaccinated with one of seven vaccines over the course of several days. The data set has been collected using microarrays. The pre-calculated results include two influenza vaccines, one containing an adjuvant and the other not. For details how to get this table, see the online tmod manual here, section “Transcriptional responses to vaccination”.

tt <- readRDS("data/vaccination_toptable.rds")

The subset contains results for 5 time points for two vaccines (10 contrasts in total). The two vaccines are denoted with “A” (no adjuvant) and “F” (adjuvant). We will run gene set enrichment on each of these time points using tmodCERNOtest.

contrasts_v <- paste0(rep(c("F", "A"), each=5), "_", 1:5)
sort_cols <- paste0("qval.", contrasts_v)
res <- map(sort_cols, ~ {
  tmodCERNOtest(tt$GENE_SYMBOL[ order(tt[[.x]]) ])
})
names(res) <- contrasts_v

The figure below shows the results.

ggPanelplot(res, filter_row_q = 1e-3, filter_row_auc = .85, q_thr=0.01)

Do these gene sets contain genes that go up or down? For this, we need to determine how many genes in each gene set are significantly going up or down.

pval <- tt %>% select(starts_with("qval"))
lfc  <- tt %>% select(starts_with("logFC"))
sgenes <- tmodDecideTests(tt$GENE_SYMBOL, lfc=lfc, pval=pval)
names(sgenes) <- contrasts_v
ggPanelplot(res, sgenes=sgenes, filter_row_q = 1e-3, filter_row_auc = .85, q_thr=0.01)

Concordance / discordance plots

Are responses on day 1 in Fluad similar to responses on day 1 in Agrippal?

  • test for interaction, or
  • use disco plots

\[d_s = \log_2 FC_A \cdot \log_2 FC_B \cdot |\log_{10}p_A + \log_{10}p_B|\]

Bibliography

Mick, Eran, Jack Kamm, Angela Oliveira Pisco, Kalani Ratnasiri, Jennifer M Babik, Gloria Castañeda, Joseph L DeRisi, et al. 2020. “Upper Airway Gene Expression Reveals Suppressed Immune Responses to SARS-CoV-2 Compared with Other Respiratory Viruses.” Nature Communications 11 (1): 1–7.
Weiner, January, David JM Lewis, Jeroen Maertzdorf, Hans-Joachim Mollenkopf, Caroline Bodinham, Kat Pizzoferro, Catherine Linley, et al. 2019. “Characterization of Potential Biomarkers of Reactogenicity of Licensed Antiviral Vaccines: Randomized Controlled Clinical Trials Conducted by the BIOVACSAFE Consortium.” Scientific Reports 9 (1): 1–14.